      subroutine RHF(S,Hcore,NORB,NATOM,Num_ele,U,E_ele)
      implicit none
      integer i,ii,j,k,l,NORB,Num_ele,NATOM,guess1
      Real*8 Hcore(NORB,NORB)       !One-electron integral 
      Real*8 S(NORB,NORB)           !Overlap matrix
      Real*8 S_DIAG(NORB,NORB)      !S^(-0.5)
      real*8 U_matrix(NORB,NORB)    !Equal to Overlap matrix
      Real*8 U(NORB,NORB,NORB,NORB) ! atom Two-electron integral
      !*****************add**********************
      character*25 file_name
      real*8 X_transfer(NORB,NORB),X_adjoint(NORB,NORB)
      real*8 P_old(NORB,NORB) ,P_new(NORB,NORB),G_ij(NORB,NORB)
      real*8 Fock(NORB,NORB),Fock_x(NORB,NORB)
      real*8 C_prime(NORB,NORB),C_new(NORB,NORB)
      real*8 E_nu,E_ele,E2,delta
      real*8 P_den(NORB,NORB),orb_energy(NORB)
      real*8 Fock_xx(NORB,NORB)
      real*8 J_martix(NORB,NORB),K_martix(NORB,NORB)
      logical logical_e
      integer count_iter
      !***************DIIS********************
      Integer,Parameter :: ndiis=7
      integer inv_DIIS,kq,kp,infodiis
      integer ipivdiis(ndiis+1)
      real*8 sumdiis
      real*8 fockdiis(NORB,NORB,ndiis),errmat(NORB,NORB,ndiis)
      real*8 cofdiis(ndiis+1),Bdiis(ndiis+1,ndiis+1)
      integer ierr,lwork
      real*8,allocatable :: work(:)
      real*8,allocatable :: ipiv(:)
      real*8 FPS(NORB,NORB),SPF(NORB,NORB)
      !***************End***********************
     !CALL RD_GEOM(NATOM,GEOM,ATOMCHG,iERRO,file_name)
     !CALL INTXC_ALLOCATE(NATOM)
     !CALL BASIS_INT(NATOM,NORB,GEOM,ATOMCHG,S,Hcore,U,baselable)
     !CALL INTXC_DEALLOCATE()
!!!==============================
        inv_DIIS=1
     ! write(*,*)'invoke RHF'
      U_matrix=S
     !do i=1,NORB
     !   do j=1,NORB
     !      do k=1,NORB
     !         do l=1,NORB
     !          WRITE(27,*)i,j,k,l,U(i,j,k,l)
     !         enddo
     !      enddo
     !   enddo
     !enddo

     !do i=1,NORB
     !       !write(*,*)Hcore(i,:)
     !       write(*,*)S(i,:)
     !enddo 
      CALL DIAG(U_matrix,S_DIAG,NORB) 
      do i=1,NORB
               S_DIAG(i,i)=S_DIAG(i,i)**(-0.5d0)
      enddo
      X_transfer = matmul(U_matrix,S_DIAG)
      X_adjoint =transpose(X_transfer)
      guess1=1 !!!Set Fock martix initial guess
      call initial_guess(guess1,NORB,X_adjoint,X_transfer,Hcore,S,P_old,Num_ele,E2)
      count_iter=0
      !P_old=0.0d0 
      logical_e=.true.
      !E2=0.0d0
      Fock_xx(:,:)=0.0d0
      !========================
      !  SCF CYCLE
      !========================
      !write(*,*)'======='
      !do ii = 1,200
      do while(logical_e)
        count_iter=count_iter+1
        !call G_matrix1(P_old,ERIall,NORB,neri,G_ij)
        call J_and_K(P_old,U,NORB,J_martix,K_martix,G_ij)
        Fock = Hcore + G_ij
        if(inv_DIIS==1)then !invoke DIIS
               ! write(21,*)'invoke DIIS'
                fockdiis(:,:,mod(count_iter,ndiis)+1)=Fock
            !   do i=1,ndiis
            !           write(100,*)Fock(i,:)
            !   enddo
            !   write(100,*)'===='
            !   do i=1,ndiis
            !           write(100,*)Fockdiis(i,:,2)
            !   enddo
               ! write(*,*)Fock(3,1),fockdiis(3,1,1) 
               ! stop
                !error matrix:ei=FPS-SPF
                !FPS
                FPS=matmul(matmul(Fock,P_old),S)
                !SPF
                SPF=matmul(matmul(S,P_old),Fock)
                !error_matrix=FPS-SPF
                errmat(:,:,mod(count_iter,ndiis)+1)=FPS-SPF
                !construct B matrix
                if(count_iter>=ndiis)then
                 !do i=1,ndiis
                 !      write(99,*)fockdiis(:,:,i)
                 !enddo
                  Bdiis=0.0d0
                  do kp=1,ndiis
                     do kq=1,kp
                        Bdiis(kp,kq)=sum(errmat(:,:,kp)*errmat(:,:,kq))
                        !if(kp/=kq)Bdiis(kq,kp)=Bdiis(kp,kq)
                     enddo
                  enddo
               !do i=1,ndiis+1
               !        write(99,*)Bdiis(i,:)
               !enddo
               !write(99,*)'======'
                !stop
                Bdiis(ndiis+1,:)=-1.0d0
                Bdiis(ndiis+1,ndiis+1)=0.0d0
                cofdiis=0.0d0
                cofdiis(ndiis+1)=-1.0d0
                !solve New matrix of coefficients
                lwork=-1
                allocate(ipiv(size(cofdiis)))
                do ii = 1,2
                        !lwork=-1
                        allocate(work(abs(lwork)))
                        call dsysv('L', ndiis+1, 1, Bdiis, ndiis+1,ipiv, cofdiis,ndiis+1, work, lwork, ierr)
                        lwork = nint(work(1))
                        deallocate(work)
                end do
                if(ierr/=0)stop 'Solve cofdiis failed!!!'
                deallocate(ipiv)
                !write(*,*)cofdiis
                !stop
              !call
              !dgesv(ndiis+1,2,Bdiis,ndiis+1,ipivdiis,Cofdiis,ndiis+1,infodiis)
              !if(infodiis/=0)stop 'Solve cofdiis failed!!!'
                !set new fock matrix
                Fock=0.0d0
                do i=1,ndiis
                   Fock=Fock+cofdiis(i)*fockdiis(:,:,i)
                enddo
                endif
        endif

       !do i=1,NORB
       !write(*,*)Fock(i,:)
       !enddo
        Fock_x = matmul(X_adjoint,matmul(Fock,X_transfer))
        
        Fock_xx=Fock_x
        !now,Fock_prime output is an eigenvector fock_prime=C'
        call DIAG(Fock_xx,C_prime,NORB)
        
        !C=X*C' here,X:X_transfer C':eigenvector
        C_new=matmul(X_transfer,Fock_xx)

        !Update the P_density matrix by new coefficient(C_new)
        call P_n_den(C_new,P_new,NORB,Num_ele)
        
        !Calculated energy
        call calc_E0(Hcore,Fock,P_new,NORB,E_ele)
        !write(*,*)E_ele
        call delta_P(NORB,P_old,P_new,delta)
        !write(28,*)delta,E_ele-E2
        if((delta < 10E-4).and.(abs(E_ele-E2) < 1.0E-6))then
        !if((abs(E_ele-E2) < 10.0E-6))then
               logical_e=.false.
              !write(21,*)'Density matrix'
              !do i=1,NORB
              !                 write(21,*)P_new(i,:)
              !enddo
              !write(21,*)'****************************'
              !write(21,*)'SCF cycle converged'
              !write(21,*)'Number of iterations:',count_iter
              !write(21,*)'Electronic energy:   ',E_ele
              !write(21,*)'Nuclear energy:      ',E_nu
              !write(21,*)'Total energy(HF):    ',E_ele+E_nu
        endif
        P_old = P_new
        E2 = E_ele
      enddo
       !write(*,*)'ok'
      endsubroutine
